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ABSTRACT 

Context. The four-planet system aroung GJ 581 has received attention because it has been claimed that there are possibly two addi- 
tional low-mass companions as well - one of them being a planet in the middle of the stellar habitable zone. 

Aims. We re-analyse the available HARPS and HIRES Doppler data in an attempt to determine the false positive rate of our Bayesian 
data analysis techniques and to count the number of Keplerian signals in the GJ 581 data. 

Methods. We apply the common Lomb-Scargle periodograms and posterior sampling techniques in the Bayesian framework to es- 
timate the number of signals in the radial velocities. We also analyse the HARPS velocities sequentially after each full observing 
period to compare the sensitivities and false positive rates of the two signal detection techniques. By relaxing the assumption that 
the radial velocity noise is white, we also demonstrate the consequences that noise correlations have on the obtained results and the 
significances of the signals. 

Results. According to our analyses, the number of Keplerian signals favoured by the publicly available HARPS and HIRES radial 
velocity data of GJ 581 is four. This result relies on the sensitivity of the Bayesian statistical analysis techniques but also depends on 
the assumed noise model. We also show that the radial velocity noise is actually not white and that this feature has to be accounted 
for when analysing radial velocities in a search for low-amplitude signals corresponding to low-mass planets. 

Conclusions. Using the Bayesian analysis techniques, we did not obtain any false positives generated by noise when analysing subsets 
of the HARPS velocities. This suggests that any signals we detect using the Bayesian detection criteria are genuine signals present in 
the data. However, some such signals could be misinterpretations caused by an insufficiently accurate noise model or by (systematic) 
noise when independent confirmation of the signals with other data sources is not possible. 

Key words. Methods: Statistical, Numerical - Techniques: Radial velocities - Stars: Individual: GJ 581 



1. Introduction 

The nature of the planetary system around GJ 581 has been dis- 
cussed in the literature by several authors and the preferred or- 
bital properties and number of planets in the system have been 
revised whenever new data has become available. While such ac- 
cumulation of knowledge is the very essence of science, contra- 
dicting results, especially those based on robust statistical argu- 
ments, have seen the system around GJ 581 become a planetary 
system with four established planets in its orbits and another two 
hypothetical ones that have been shown to exist in some studies 
but have been shown detections of false positives in others. 

Astronomers first started counting the planets orbiting GJ 
581 in 2005 when Bonfils et al. (2005) reported they had found 
a Neptune-mass planet (m p sini = 16.6M ffl ) with an orbital pe- 
riod of roughly 5.4 days. This discovery, made using the High 
Accur acy Radial Velocit y Planet Searcher (HARPS) spectro- 
graph dMavor et all 120031) mounted on the 3.6 m telescope at 
the European Southern Observatory (ESO) at La Silla, Chile, 
was soon accompanied by discoveries of two additi onal super- 
Earth -mass planets made by the same team in 2007 dUdrv et all 
120071) . These two additional planets were reported to have min- 
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imum mas ses of 5.0 and 7 .7 M© and orbital periods of 12.9 and 
84 davs. lUdrvetal.1 d2007l) argued that the outer planet (GJ 581 
d) was on a slightly eccentric orbit (e =0.20+0. 10) but also noted 
that its orbit could be consistent with a circular one as was the 
case for the other two planets. 

In 2009. lMavor et all d2009h reported that additional HARPS 
velocities revealed the existence of a fourt h planet orbi ti ng GJ 
581 with and orbital period of 3.1 days. Ma yor et alj ([2009) 
detected this planetary companion from 119 velocity measure- 
ments and claimed it to have a low minimum mass of 1.9 M®. 
They also revised the orbital period of GJ 58 1 d to 66.4 days and 
stated that the period of 84 days obtained bv lUdrv et all (120071) 
was in fact a yearly alias of the 66.4-day signal. This is a text- 
book example of problems aliasing can cause to detections of 
periodic signals in radial velocities and other time series. 

Counting the planets orbiting G J 581 got complicated 
in 2010 and 2011 when IVogtet al.l d2010h claimed that a 
combined d a taset with the 119 HARPS measurements of 
iMavor et al.l ([2009) and anot her set of 122 H igh Resolution 
Echelle Spectrograph (HIRES: IVogt et all[T994l) velocities from 
the 10 m Keck-II telescope at Hawaii actually contained the 
periodic signatures of - not only the four planets observed by 
Mayor et al. (2009) - but also two other signals corresponding to 
low-mass planets with minimum masses of 7.0 and 3.1 M ffi and 
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orbital periods of 433 and 36.6 day s, respectively. The immedi- 
ate conclusion of IVogt et al.l (12010 1) was that the latter might be 
a good candidate for a habitable planet with its Earth-like mass 
and its orbit in the middl e of the stellar hab itable zone of G J 581. 
However, the analyses o flTuomil d2011l) and lGregorvl(l201 ll) soon 
revealed that the existence of this candidate habitable world was 
not supported by the data and that the evidence in favour of the 
outer pl anet with an orb ital period of 433 d ays was at most sug- 
gestive. Gregory] d201 lb and ITuomil d20 lib concluded that there 
were five and four planets, respectively, in the system depending 
on the subjective choices of the probabilistic detection threshold 
and the prior densities of the parameters of the statis tical mod- 
els. Ye t, the existence of the four planets detected bv lMavor et al.l 
(2009) was easily verified using th e Bayesian statistical tests of 
ITuomil d2011l) and lGregorvl d2011l) . 

Other author s have also r e-analy sed t he HARPS and/or 
HIRES data of iMavor et al.l d2009l) and IVogt et al.l d2010h 
and received evidence that supports the existence of 
the hypothetical habitable zone pl anet GJ 581 g (e.g . 
Angl ada-Escu de & Dawsonl 120111 : iTadeu dos Santos et al.L 
2012). However, it has been clear from the start, as also stated 
by almost every author discussing the system and estimating 
the number of planets in it, that only additional data can firmly 
verify or falsify the existence of the two proposed companions 
GJ581fandg. 

Such data was made available bv lForveille et al.1 (1201 ll) and 
increased the amount of HARPS data of GJ 581 by another 121 
velocities. This increase made the total number of HARPS preci- 
sion velocities 240. However, despite attempts of recovering the 
proposed planets f and g from this data using sta ndard Lomb- 
Scargle periodograms ( Lomb. 1 1 976c IScargleL 1 1 9821) and^ 2 min- 
imisations, Forveille et al.l (1201 lli concluded that these two plan- 
ets could not exist orbit ing the star wit h the masses and orbital 
properties reported by dVogt et alll2010 ). Recently, a re-analysis 
of these HARPS velocities backed up by dynamical arguments 
was found to support the existence of a habitable-zone super- 
Earth with an orbital period of 32 days and a minimum mass 
of 2.2 M ffi JVogt et ul.. 2012) but another study taking into ac- 
count the noise correlations did not find evidence in favour of 
the existence of this signal in the same data and also casted 
doub t on the sig nificance of the signal corresponding to GJ 581 
d(|Baluevl[2012). Therefore, considering the different estimates 
for the number of planets (parameter k) favoured by the data in 
the very recent history suggests that this debate is far from over. 

If of planetary origin, any periodic variations in radial veloci- 
ties should be present in data obtained using different telescope- 
instrument combinations. The availability of such independent 
confirmation is essential to exclude spurious signals that can be 
caused by instrument instability and insufficient noise modelling 
of the features in the measurement noise of the instrument. 

Using the combined set of HARPS and HIRES velocities of 
GJ 581, we perform a global search for low-a mplitude signals 
in an attempt to see whether the conclusions of Forveill e et al.l 
d20"Tll) . stating that k — 4, hold for results obtained us- 
ing Bayesian analysis techniques (e.g. ITuomil 1201 ll [2012; 
Tuomi et al], 1201 It iTuomi & JonesL 120121). We also re-analyse 
the 240 HARPS velocities of IForveille etall (1201 ll) to (1) en - 
able comparison o f our results wi t h thos e ofJBonfils et al. (2005), 
Udrv et al l (120071) . IMavor et al.1 (120091) . IForveille etaiTlhoilh . 
Vogt et alT(l2012l) . and Baluevl d2011l) ; to (2) quantify the sen- 
sitivity of the Bayesian detection methods w.r.t. the frequen- 
tist periodogram tools used throughout the literature; and to (3) 
see whether the apparent sensitivity of the Bayesian detection 
thresholds is prone to false positives in practice. Our goal is to re- 



ceive a reliable estimate for k in the case of GJ 58 1 radial veloci- 
ties and to determine the significances and properties of possible 
additional signals if k > 4. In our Bayesian analyses, we use the 
"standard" white noise models assuming independent and iden- 
tically distributed measurements with Gaussian white noise but 
we also test models with red noise features, i.e. we model the 
data using a moving average (MA) model that accounts for the 
noise correlations i n the data in a time- scale of few dozen days 
(e.g. lBaluevLl20H ITuomi et al.L l2012allbh . 

2. Statistical analysis 

2.1. Bayesian statistics 

When assessing the optimal value of k, i.e. the number of signals 
statistically significantly present in a time series, we take advan- 
tage of Bayesian model selection that enables the simultaneous 
comparison of any number of models - in this particular case, 
models with different k. We denote these models as Mk, for all 
k = 0, . ..,p, and calculate the corresponding probabilities of each 
model in a standard way as 



P(M k \m) 



P{m\M k )P(M k ) 



(1) 



where m denotes the measurements; P(m\Md is the Bayesian 
evidence, or probability of obtaining the data if the model Mi 
indeed was a "correct" model; and P{Mi) is the prior probability 
of the model At, that quantifies the subjective level of confidence 
we have on the validity of this model a priori to analysing the 
measurements. As can be seen in Eq. (G3, the probabilities are 
scaled in such a way that the sum of the individual probabilities 
of models Mk, for k = 0, p, is unity indicating that any poste- 
rior probabilities obtained using this equation are dependent on 
the selected set of candidate models and only reflect the relative 
probabilities of the models in this set. 

Finding the optimal value of k based on model selections us- 
ing Eq. dTJ, however, is not as such a sufficient way of counting 
the periodic signals in a given data set. While a given k might 
receive the greatest posterior probability, the claim that there re- 
ally are k signals in the data should not be made lightly. This is 
especially the case if the probability of a model with some lower 
value of k receives estimates of similar magnitude. For this rea- 
son, we require that the probability of a model Mk+\ has to be 
at least 150 times that of a model Mk to conclude that there is 
evidence in favour of the existence of k + 1 Keplerian signals 
in the data. This threshold corresponds to an inte rpretation of 
"strong evidence" given bv lKass & Raftervl (1 19951) to make de- 
cisions based on data. While completely subjective, we adopt 
this threshold but do not cons ider signals t o exist in the data un- 
less the other two criteria of ITuomil d2012h are also satisfied, i.e. 
that the estimates of the velocity amplitude parameter K are sta- 
tistically significantly greater than zero and that the period P is 
constrained from above and below. If all these criteria are satis- 
fied, we state that the existence of a signal is supported by data. 

While the prior probabilities in Eq. dTJ are only simple num- 
bers, the Bayesian evidences are in fact complicated integrals 
over the whole parameter space of model parameters that de- 
pend on the value of the product of a likelihood function and a 
prior probability density of the model parameters, i.e. the formu- 
lation of the statistical model and the physical interpretation of 
its parameters given the available data. For a given model, this 
integral, sometimes called the marginal integral, is of essence in 
any model comparison problems in the Bayesian context and the 
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ability to compare models in a probabilistic way depends on the 
ability to estimate this integral. We write this integral as 



P(m\M k ) 



-j 



l{m\e k ,M k )n(6 k \Mk)dO k , 



(2) 



where the parameters of the kth model, 6 k , in the parameter space 
Q. k , are interpreted using the likelihood function / and their prior 
probability density n that describe the mathematical details as 
well as the physical interpretation of the statistical model M k . 
This is also the reason we write the likelihood and the prior as 
conditioned on the model. 

We estimated the integral in Eq. (O usin g the truncate d 
posterior mixture (TPM) estimate of Tuomi & Jones! d2012l) . 
This estimation requires the availability of a statistically rep- 
resentative sample from the parameter posterior density and 
we drew such samp les using the adaptive Metropolis algorithm 
dHaario et alll200lh that is simply an adaptive version of the fa- 
mous Metr opolis-Hastings Markov chain Monte Car lo (MCMC) 
algorithm dMetropolis et all 119531: lHastingsl 1 1 9701). This sam - 
pling technique has been used succesfully in e.g. Tuomil (l2012h . 

Apart from the chosen set of candidate models (i.e. the def- 
inition of their likelihood functions), the only other subjective 
issue in Bayesian statistics is the prior information quantified 
using prior probability densities of the model parameters and 
prior probabilities of the different models. Standard (frequentist) 
analyses of radial velocity data actually incorporate prior infor- 
mation as well, which makes them somewhat Bayesian in reality, 
but do not easily enable testing different prior opinions in prac- 
tice. For instance, when searching for a solution using typical^ 2 
minimisations, the amount of excess noise in the measurements, 
or parameter cry, is commonly fixed to an a priori selected value 
(e.g. ctq), which corresponds to setting the prior probability of 
this parameter equal to a delta-function as ji{ctj) = 6(<tj - o~o). 
This prior assumes that only a value of cry = ctq is a possible 
one and ignores the rather likely possibility that this parame- 
ter might not be equal to <tq. If this prior choice was made in 
the Bayesian context, any scientifically oriented reader would 
immediately dismiss this choice as too limiting and would not 
consider the obtained results trustworthy. However, such prior 
limitations, i.e. setting prior densities to zero in chosen subsets 
of the parameter space, are commonly used in statistical anal- 
yses but should be used with care and only when the physical 
interpretation of the parameters justifies the limitations. 

Similar "hidden" prior assumptions are made throughout the 
literature announcing detections of exoplanets around nearby 
stars. For instance, finding a solution to some data by assum- 
ing that planetary eccentricities are equal to zero represents a 
similar prior choice. It could be argued that in "dynamically 
packed" systems only eccentricities close to zero are viable by 
enabling the long-term stabilities of these systems, but the pos- 
sibility that some non-zero orbital eccentricities could also pro- 
vide stable systems cannot be ignored because enabling eccen- 
tricities as free parameters effectively increases the number of 
free parameter in the model by 2k compared to a model with 
fixed eccentricities and this can have considerable effects on the 
significances of the obtained solutions. As fully Bayesian, these 
hypotheses could be easily tested by treating the model with ec- 
centricities fixed to zero as another model in the set of candidate 
models. In practice, this can be seen as a comparison of different 
prior densities of the model parameters. 

The argument of hidden prior information can in fact be 
taken even further. Any non-Bayesian analyses correspond to 
Bayesian ones with the assumption that the prior probability 
densities of model parameters are actually "flat", i.e. uniform 



distributions in the parameter space. This is a typical uninfor- 
mative choice for a prior density and is applied widely in the 
statistical literature to various parameter estimation problems. 
However, it still is a subjective choice and it does have an effect 
on the obtained results. For instance, a typical non-linear trans- 
formation of the parameter vector from 9 to ff reveals the conse- 
quences of this prior choice. If the prior for parameter 6 is a uni- 
form density, that for ff has to be chosen by using the Jacobian 
of the corresponding transformation to receive the same results. 
If the prior of ff is chosen uniform as well, the obtained anal- 
ysis results would not be the same for these two parameterisa- 
tions but would depend heavily on the chosen transformation, 
which emphasises the fact that prior densities are important en- 
tities that cannot be ignored when performing statistical analy- 
ses. This in turn implies that the only logically consistent frame- 
work of statistics is the Bayesian one. While these issues are 
discussed widely in the statistical literature, we refer to the ex- 
cellent book Sta tistical Decisi on Theory and Bayesian Analysis 
bv .1. O. Berger dBergerl[l980h . 

Because priors are subjective, we present our results only 
based on the priors we choose. In practice, we adopt the same 
pr ior probabili ty densities of model parameters that were used 
in ITuomil d20 1 2h . We also use the same prior model probabilities 
as ITuomil d2012l) to take into account the fact that we believe 
finding k+ 1 planets in any given system is less likely than finding 
k of them. Because posterior samplings can be used to obtain 
estimates of the parameter probability densities that enables the 
computation of any point- and uncertainty estimates, we report 
all the results using the maximum a posteriori (MAP) estimates 
and the correspond ing 99% Bayesian c redibility sets (BCSs; as 
defined in e.g. lTuomi & Kotirantall2009l) . 



2.2. Periodograms 

Because it is useful to calculate the periodograms of time se- 
ries, especially when searching for periodic signals, we discuss 
the resulting power spectra as well for the GJ 581 velocities. 
Identifying the strongest powers in the periodogram , such as 
the fa mous Lomb-Scargle periodogram (Lomb, Il976t IScargle. 
1 19821) . helps choosing the initial states of posterior samplings by 
pointing towards the regions in the period space where periodic 
signals are the most likely to be found. The ability to use infor- 
mation from periodogram analyses can help decreasing the burn- 
in period, i.e. the length of the Markov chain before it has con- 
verged to the global solution, of MCMC samplings considerably. 
However, because periodogram analyses are commonly used to 
assess the number k in practice, we first discuss the caveats of 

such assessments; 

As noted bv ITuomil d2012l) . periodograms and related anal- 
ysis methods rely on the assumption that the model residuals 
have been calculated correctly. For instance, when searching for 
a k+ 1th signal in a noisy time series, the residuals are commonly 
calculated by assuming that there are in fact only k periodicities 
in the data. This obvious contradiction of the assumption and the 
statistical test performed to search for another signal is not only 
made throughout the literature when searching for planets using 
the radial velocity technique, but its effects on the detectability 
of periodic signals have not received as much attention in the 
same literature as they should. A clear example is the case of 
HD 10180 that was reported to hos t six, possibly seven, plan- 
ets orbiting it by lLovis et ail (1201 ll) . They used power spectra 
in their attempt of assessing the number k for this star and fell 
a victim of the above contradiction of the statistical test and the 
assumption this test is based on. Contrary, the Bayesian analysis 
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of the same exact data performed by Tuomj| (1201 2l) concluded 
that the preferred estimate for k is actually as high as nine. This 
example only demonstrates the weakness of the periodogram- 
based methods in assessing k but can also be seen as an obvious 
consequence of the assumption: if it is assumed that there are k 
signals in a data set, the statistical tests performed based on this 
assumption will certainly be biased in favour of this assumption 
being true and prevent the detection of a k + 1th signal unless it 
is a strong one and clearly present in the data. 

In an attempt t o overcome the above shor tcomings of peri- 
odogram analyses, Anglada- Escude & Tuomil d2012l) proposed 
a gene ralised version of the periodogram defined by ICum ming 
(120041) for the purpose of detecting multiple periodicities in time 
series . This recursive periodogram (Anglada-Escud e & Tuomil 
[2012b can be used to obtain the power spectrum for model resid- 
uals by simultaneously adjusting the parameters of the previ- 
ously detected k signals for each test frequency of the k + 1th 
one. Clearly this method is not prone to such a severe contradic- 
tion as the traditional periodograms but it is still not free of its 
own conseptual problems such as the fact that one has to move 
sequentially by increasing k one step at a time when searching 
for several periodicities. While certainly an improvement to the 
classical periodograms, the consequences of hidden prior infor- 
mation cannot be addressed using this modified periodogram ei- 
ther. 

We report the results of standard Lomb-Scargle periodogram 
analyses together with those from the Bayesian analyses. We 
also report the resulting false alarm probabilities (FAPs) of the 
signals we detect in the GJ 581 velocities obtained using the 
HARPS spectrograph. 



3. Sequential analysis of HARPS velocities 

The HARPS velocities provide an interesting opportunity to 
compare the sensitivity of the Bayesian analysis methods to the 
traditional periodogram ones used throughout the literature. For 
this reason, we construct a timeline for GJ 58 1 where we analyse 
the available HARPS data afte r each yearly obser v ing period. 
We als o an alyse the data sets of Bonfils et al.1 (l2005l) . IUdrv et al.1 
(2007), and M ayor et al.l(l2009l) and compare the best estimate of 
k with an estimate that could have been received had Bayesian 
tools been applied to these velocities since 2005. 

We show the number of HARPS velocities after each ob- 
serving period and at the time when the three studies describ- 
ing detections of the four planets orbiting GJ 581 and reporting 
the data were published (these three stud i es are : iBonfils et all 
12005b lUdrv eTall l2007t iMavor & Ouelozl Il995l) . The numbers 
of measurements and data baselines of these subsets of HARPS 
data are shown in Table Q] In the following subsections we de- 
scribe briefly the results of our analyses of the subsets of HARPS 
data in Table Q] We note, that we assign numbers to the signals 
we detect by calling the shortest periodicity as signal 1, the sec- 
ond shortest as signal 2 and so forth. We also refer to the signals 
corresponding to th e four planet cand idates by using the letters 
b, c, d, and e as in lMavor et all d2009l) . 



3.1. Bonfils etal. data 



Table 1. Number of signals detected from the subsets of HARPS 
data as a function of time as reported in the literature (ki), using 
Lomb-Scargle periodograms (£2), and using Bayesian tools (£3). 
The corresponding data baselines (T b s ) and numbers of data (AO 
are also shown. T he subsets correspond to the published data 
sets and subsets of Forveille et al. data set after each full 

observing period (FOP). 



Data subset 


N 


r obs [days] 


h 


k 2 


h 


Bonfils etal. (2005) 


20 


440 


1 


1 


1 


First FOP 


25 


457 




2 


2 


Second FOP 


43 


827 




2 


2 


Udrv et al. (20071 


50 


1050 


3 


3 


3 


Third FOP 


76 


1197 




4 


4 


Mayor et al. (2009) 


119 


1570 


4 


4 


4 


Fifth FOP 


131 


1904 




4 


4 


Sixth FOP 


197 


2312 




4 


5 


Forveille etal. (2011) 


240 


2543 


4 


4 


5 



When analysing the IBonfils et al.l d2005l) data set, our posterior 
samplings could rapidly verify that it contains the clear plane- 
tary signature of GJ 581 b. This signal was seen easily with pos- 
terior samplings (with P(k = 1) g 5.9 x lO lo P(k = )) and pe- 
riodograms (FAP < 0.0001) and IBonfils etaf] d2005l) were also 



very confident about its existence based on their analyses that 
we assum43 relied on periodograms and^- 2 minimisations. 

Our samplings of a two-Keplerian model also spotted hints 
of a second periodic signal between roughly 10-16 days. The 
two-Keplerian model received a significantly greater posterior 
probability (by a factor of 640) than the one-Keplerian one, 
though the samplings could not constrain the period of the sec- 
ond signal very accurately. Also, parameter K of this second 
signal did not differ statistically from zero, which forced us to 
conclude that there i s only one signal significantly present in the 
BonfilsetalJ (120051) data set. Yet, this signal exists with little 
doubt because taking it into account increases the model proba- 
bility by a factor of l.lxlO 11 compared to the model with k = 
despite the low number of measurements (see TableQ}. 



3.2. First full observing period 

The end of the first full observing period provided only five 
more spe ctra and cor r espon ding velocities of GJ 581 than the 
data of (Bonfi ls et all |2005). However, this changed the situa- 
tion with respect to the second signal considerably. The two- 
Keplerian model became 4.6xl0 4 times more probable than the 
one-Keplerian one, and parameter K of the second signal was 
significantly greater than zero with K2 = 3.6 [1.5, 5.7] ms _1 . 
Despite the fact that the period of this signal had a density with 
several modes between roughly 12 and 16 days, it was well con- 
strained from above and below and thus satisfied all the detec- 
tion criteria. The second periodic signal was also detected by 
the periodograms and the corresponding FAP of this signal was 
0.0075. 

The probability densities of the velocity amplitude and pe- 
riod of the second signal are plotted in Fig.Q] The multimodality 
of the period is clear from the top panel in Fig. Q] but the period 
is still constrained from above and below and its 99% BCS is 
[11.98, 16.27] days. We note that the MAP estimate of the pe- 
riod is 12.9 days. 

We sampled a three-Keplerian model as well, but unsurpris- 
ingly a third periodic signal did not get constrained by the avail- 
able 25 velocities. 



1 In fact, Bonfils et al. ( 2005) do not mention the statistical methods 
they used to obtain their solution. 
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Fig. 1. Distributions of the velocity amplitude and period of the 
second signal (corresponding to GJ 581 c) as obtained from the 
partial HARPS data after the first full observing period. The 
solid curve shows a Gaussian function with the same mean and 



- Mode: 3.1506436 : 3-7079022 

W: 3.1503417 //; 18.938814 
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1.0216014 //: 0.17077295 

0.31918028 




Fig. 2. As in Fig. Q] but for the velocity amplitude and period of 
the signal with the shortest pe riodicity (corresponding to GJ 581 
e) from the lUdrv et all (120071) data. 



variance. 



3.3. Second full observing period 

At the end of the second full observing period, the amount of 
velocities had increased to 43, and we could easily sample mod- 
els with k = 0, 3. For this data set, the three-Keplerian model 
implied that there are at least three periodicities in the data but 
the period of the third signal turned out to have a density with 
three modes. These modes were found at periods of 42.3, 67, and 
84 days. We could not set constraints to the period, and could 
not be sure whether our Markov chains had indeed converged 
to this multimodal density - but assuming they had, the three- 
Keplerian model would have had a probability of 620 times that 
of the two-Keplerian one whose parameter space was well sam- 
pled and provided two clear periodicities at 5.37 and 12.9 days. 

When we sampled the parameter space of a four-Keplerian 
model, however, we could identify a solution for a fourth pe- 
riod at 3.1 days. This solution, while only suggestive because 
the period was not constrained from above or below but was only 
apparent as a global maximum in the otherwise "flat" probabil- 
ity density, actually helped constraining the third periodic signal 
from above and below according to our detection criteria. While 
this subset of HARPS data could not be used to claim a detec- 
tion of four periodicities, the results suggested that only a modest 
increase in the number of available data would be likely to pro- 
vide detections of four signals. Yet, according to our detection 
criteria, we could spot only two signals confidently in this data 
set, namely, those with periods of 5.37 and 12.9 days. These two 
signals were also detected by the periodogram analyses (Table 

ID- 

3.4. Udry et al. data 

lUdrv et al.1 (120071) reported that the iBonfils et all (120051) data 
already actually suggested there might be an additional com- 
panion orbiting GJ 581 with an orbital period of 13 days with 
a modest significance. With additional 30 data points and im- 
proved barycentric corrections, their periodogram analyses and 
X 1 minimisations pointed towards two new periodic signals, cor- 
responding to two new planets orbiting GJ 581 with orbital peri- 
ods of 13 and 84 days. 

These two s ignals were indeed significantly present in the 
lUdrv et al.1 (|2007) data according to our Bayesian analyses. 
However, we could go further than that and identify a fourth 
periodicity at 3.14 days corresponding to the fourth planet can- 
didate reported bv lMavor et al.l (1200 9). We show the probability 
densities of velocity amplitude and period of the fourth signal 
(Fig. 13 to demonstrate that the fourth periodicity was detected 
as a clear probability maximum but it did not quite satisfy the de- 
tection criteria, namely, parameter K of the 3.14 day signal was 



not significantly greater than zero by having a tail that extended 
all the way to zero in the amplit ude space (Fig. [2] top panel). 
However, the signals reported bv lUdrv et alj d2007l) were clearly 
detected in the data. 

The periodogram analyses revealed the existence of three 
signals in this data set as well. The third most significant signal 
was found at 82.6 days with a FAP of 0.001 but the interpretation 
of the origin of this signal corresponding to the repo r ted p eriod 
of 84 days was less clear than claimed bv lUdrv et alj d2007l) be- 
cause the sixth stronges power in the window function was sus- 
piciously close to this period at 87 days making it possible that 
this periodicity could be due to sampling and not to a genuine 
signal. 

3.5. Third full observing period 

After the third full observing period with 76 radial velocities, the 
3.14 day signal that was close to satisfying the detection criteria 
in the lUdrv et al.l d2007h data set was well constrained and satis- 
fied all the detection criteria. With these velocities, the period of 
the outer planet candidate received a glob al maximum at 67 days 
instead of the periodicity at 84 days in the lUdrv et al.l 12007) data 
set . Essentially, we o btained the same solution that was reported 
by iMavor et al.l d2009l) after four full observing periods and 43 
more velocities. 

To demonstrate that the solution indeed satisfied the detec- 
tion criteria, we show the posterior densities of the period and 
velocity amplitude of the "weakest" 3.14 day signal in Fig.[3](top 
panels) together with the phase-folded signal (bottom panel). 
Because the corresponding four-Keplerian solution also had a 
posterior probability that was 5.4xl0 7 times greater than that of 
the three-Keplerian model, we could conclude that the 76 epochs 
of HARPS data obtained after three full observing periods were 
already clearly in favour of a four-Keplerian model. We note that 
despite the eccentricity prior that penalises high eccentricities, 
the eccentricity of the outer companion had a MAP estimate of 
0.30, yet, its uncertainties were considerable and it was found 
to be statistically indistinguishable from zero with a 99% BCS 
of [0, 0.58]. The four signals were also obtained by the peri- 
odogram analyses and the FAP of the fourth signal was 0.001, 
which indicates that this signal was significantly present in the 
data. 

We sampled the five-Keplerian solution and detected two 
probability maxima for the fifth signal at 25 and 42 days. 
However, neither of these was found to correspond to a solution 
that satisfied the detection criteria of another signal - the corre- 
sponding periods did not get well constrained and the velocity 
amplitudes received values were not found distinguishable from 
zero. 
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Fig. 3. Probability densities of parameters P and K of GJ 581 
e as in Fig. Q] and the phase-folded MAP signal (bottom panel) 
with the other three signals subtracted. 



3.6. Mayor etal. data 

iMavor et al.l (|2009) reported the detection of the fourth and in- 
nermost planet candidate orbiting GJ 581 with an orbital period 
of 3.14 days after having observed 119 spectra of the star. This 
new data set included velocities based on all the existing spec- 
tra after improved spectral reductions and corrections account- 
ing for e.g. a changing inter nal pressure of the Th-Ar calibra- 
tion lamps caused by aging. Mavor et al. (2009) stated that they 
discovered a fourth periodic signal in the velocities using "non- 
linear minimisations with genetic algorithms" and periodogram 
analysisQ 

We c ould obtain the sam e four-Keplerian solution easily us- 
ing the IMavor et al.l (120091) velocities. However, we also tested 
whether we could spot a fifth signal in this set of 119 radial 
velocities. Indeed, sampling the parameter space of the five- 
Keplerian model identified three probability maxima for a fifth 
periodicity in the data. These maxima were found at periods of 
17, 34, and ~400 days but the corresponding Markov chains 
did not actually converge to a posterior density in the sense that 
their means and variances did not show converge as a function 
of the chain length. While this might be a result of insufficient 
chain length, we drew samples of up to few 10 s members from 
the posterior density, which should be a sufficient length given 
the acceptance rate of roughly 10% throughout the samplings 
dTuomil l2012h . Also, for the local maxima at 17 and 34 days, 
the velocity amplitude parameter was not found statistically dif- 
ferent from zero, which violates the detection criteria. Similarly, 
the global maximum of the fifth periodicity at ~400 days (like ly 
corresponding to the periodicity reported bv lVogt et alj d2010t) at 
433 days) was found to have an amplitude consistent with zero 
despite the fact that its MAP estimate was found to be as high 
as 1.7 ms _1 . The signal corresponding to this global maximum 
was also found to have an extremely poor phase-coverage due to 
annual gaps in the HARPS velocities. 



3. 7. Fifth full observing period 

The 131 velocities obtained after the fifth full observing period 
did not change the soluti on much with resp ect to the one ob- 
tained usin g the data set oflMayor et alJ (120091) . The four signals 
reported bv lMavor et al.l d2009l) were very clearly present in the 
data and we moved on to sample the parameter space of a five- 
Keplerian model. This time, we spotted four different probabil- 
ity maxima for the period of the fifth signal at 32, 42, 190, and 
420 days. However, none of these maxima corresponded to solu- 
tions satisfying the detection criteria by having either posterior 
probabilities below the detection threshold of by having velocity 
amplitudes statistically indistinguishable from zero. Therefore, 
we could not claim that the data supported the existence of more 
than four signals according to our detection criteria. 

Curiously, the 42 day probability m aximum corresponds to 
the signal reported bv lVogt et al.l (l2010 l ) whereas t he 32 day one 
corresponds to the signal reported by IVogt et al.l d2012l) . These 
two probability maxima, despite the fact that they do not sat- 
isfy the detection criteria, are therefore existing features of the 
parameter density of the five-Keplerian model. Based on these 
131 velocities, however, it is impossible to state whether these 
probability maxima correspon d to actual signals in the data or 
false p ositives in the analyses of lVogt et"al] d20 1 Ot) and Vogt et al. 
(2012). 



3.8. Sixth full observing period 

After the sixth full observing period, the total amount of HARPS 
velocities for GJ 581 had increased to 197. These data further 
im proved the constra ins of the orbits of the four planets detected 
bv lMavor et al.l d2009l) and increased their significances w.r.t. the 
131 velocities that had been available a year earlier after the fifth 
full observing period. When sampling the five-Keplerian model, 
we could identify an interesting solution for the fifth signal. 

According to our samplings, there were actually two proba- 
bility maxima in the five-Keplerian model corresponding to 170 
and 190 day periodicities for the fifth signal. These maxima, by 
being rather close to one another in the period space, enabled 
drawing a sample from the posterior density and our Markov 
chains were found to "bounce" between the two maxima yield- 
ing a sample from this bimodal solution. This solution satis- 
fied all the detection criteria by having a well constrained pe- 
riod, velocity amplitude that was significantly greater than zero, 
and a model probability that was 4.7 xlO 3 times greater than the 
probability of the four-Keplerian model, which indicates that the 
favoured value of A: is 5 for this data set. While the bimodality 
is a rather unsatisfactory feature of a probability density - es- 
pecially when the goal is to quantify the properties of signals 
corresponding to candidate planets - it is by no means surpris- 
ing considering that we are searching for periodic signals in the 
first place. Similar bimodal densities ha ve been obser ved before 
when dealing with radial velocities (e.g. lTuormll2012l) . We show 
the posterior density of the period of the fifth signal in Fig. [4] 
We could not find any other significant periodicities in this data. 
Also, we could not detect the fifth signal using the periodogram 
analyses of the data. 



3.9. Forveille etal. data 



2 IMavor et alj (2009) did not actually specify which measure of 
model goodness they minimised, which makes it impossible to repli- 
cate their results. However, the term "minimisations" that they used 
suggests that the chosen measure was x 1 statistics which corresponds 
to a delta-function prior for the jitter. 



iForveille et al.l (1201 ll) reported that they could spot four and only 
four periodic signals in their set of 240 HARPS velocities. This 
was also a simple task for our periodogram analyses as well as 
for the Bayesian posterior samplings and we continued to sam- 
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Fig. 4. Probability density of the fifth periodicity in the HARPS 
data after the sixth full observing period. 

pie the parameter space of the five-Keplerian model. These sam- 
plings indicated that the 190 day periodicity observed in the data 
after the sixth observing period gained greater significance while 
the 170 day period ceased to correspond to a probability maxi- 
mum. The probability of the five-Keplerian model was found to 
be 5.7xl0 4 times more probable than the four-Keplerian one, 
which shows that the significance of the fifth periodic signal at 
190 days received roughly ten times greater significance than it 
did for the data after the sixth full observing period. This peri- 
odicity was also well constrained and satisfied all our detection 
criteria. The periodogram analyses could not be used to detect a 

fifth periodogram in the data co nfidently. 

However, as also reported by Vogt et al.l d2012l) . we also spot- 
ted another probability maximum corresponding to a periodicity 
at 32 days. Treating the 32 and 190 day solutions as separate 
models, we calculated their respective probabilities. These prob- 
abilities were found to be 7% and 93%, respectively, which indi- 
cates that the global solution of the five-Keplerian model corre- 
sponds to a fifth periodicity at 190 days. Because of the existence 
of the probability maximum at 32 days, we performed global 
samplings of a six-Keplerian model but failed to find a signal 
satisfying the detection criteria. We did, however, spot a proba- 
bility maximum for the period of the sixth signal at 32 days, but 
it could not be constrained from above or below, which violates 
one of our detection criteria. We show the phase-folded fifth sig- 
nal in Fig. [5] (bottom panel). While the phase-folded signal is 
not visually very impressive, we also plotted the distributions 
estimating the period and amplitude of this signal in Fig. (top 
panels) which indicate that this signal complies with our detec- 
tion criteria. 



4. Analysis of combined radial velocities 

In this section we perform a robust analysis of the combined 
HARPS and HIRES velocities of GJ 581. The base li ne of the 
available 240 HARPS velocities of iForveille et all (1201 lb is 
2543 days. We plotted the four-Keplerian model residuals of the 
HARPS velocities in Fig.|6](top panel). Clearly, these velocities 
show the remarkable stability (with an RMS of only 1.98 ms" 1 ) 
of the HARPS spectrograph and indicate that, in addition to ob- 
vious yearly gaps, the HARPS data can be expected to provide 
accurate constraints to possible additional periodic signals in the 
data. 

The residuals of t he four-K e pleria n model of the HIRES pre- 
cision velocities of Vog t et ail (1201 Ol) are also shown in Fig. [6] 
(bottom panel). These velocities have an RMS of 2.82 ms -1 , 
which implies an impressive precision, but falls short of the 
precision of the HARPS velocities for GJ 581. The baseline of 
the HIRES data is 4000 days, which is almost twice that of the 
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Fig. 5 . As in Fig. |3]but for the fifth signal in the iForveille et al.l 
(2011) data. 
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Fig. 6. HARPS (top) and HIRES (bottom) residuals of the four- 
Keplerian model. 



HARPS one. Therefore, combining the HIRES and HARPS ve- 
locities would certainly provide greater accuracy than either one 
of them alone by having a greater overall precision and longer 
baseline. Also, the differences in the time-distribution of these 
two sets can certainly help to distinguish genuine Doppler peri- 
odicities and their aliases from one another and from spurious 
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signals generated by noise and insufficient modelling by provid- 
ing a better sampling than the individual sets. 

4. 7 . White noise model 

As we already analysed the HARPS velocities individually in 
Section 13.91 and because the HIRES velocities have been anal- 
ysed individually in e.g. Greg ory! (1201 ll) . where it was concluded 
that they contain two significant periodicities (the two strongest 
signals at 5.4 and 12.9 days), we move on directly to analyse the 
combined data set with the pure white noise model. 

Using the Gaussian white noise model, the four signals re- 
ported by Ma yor et al.l ((2009) were indeed found clearly in the 
combined HARPS and HIRES velocities and our Markov chains 
converged to the corresponding four-Keplerian solution with 
ease. All these four signals were well constrained and satisfied 
our detection criteria. 

When sampling the parameter space of the five-Keplerian 
model, we searched for signals in the period space between 12.9 
and 66.7 days and the space beyond 66.7 days up to a cutoff 
period of r max = 10r o b s , where r o b s is the data baseline, sepa- 
rately and treated them as different models. Effectively, this cor- 
responds to comparing two different prior models, i.e. in this 
case prior limitations (see Section 12.11 ), with period ranges as 
P € [12.9,66.7] and P e [66.7, r max ], where P and r max have 
the unit of days, with each other. The cutoff period exceeds the 
data baseline, but we chose such a long upper limit for the pa- 
rameter space because signals with periods excee ding the data 
baseli ne can be detected in radial velocity data (Tuomi et al.L 
2009) and this cutoff choice also enables the detection of sig- 
nals that are only present in the data as trends with no or little 
curvature. We chose the two other cutoffs of the period space 
because we expect a priori the inner system of GJ 581 to be 
"dynamically packed" in the sense that because of the existence 
of planets with orbital periods of 3.1, 5.4, and 12.9 days, it is 
likely that there are no stable orbits for additional planets with 
orbital periods less than 12.9 days. This leaves two subspaces 
of the period space where stable orbits could exist, namely, the 
ones we define above, i.e. orbits between GJ 581 c and d and 
orbits beyond GJ 581 d. 

Because of the division of the period space into two parts, 
we could draw samples from both subspaces relatively easily. 
Samplings of the parameter space of a five-Keplerian model re- 
vealed a clear probability maximum for the period of the fifth 
signal at 43.7 days. We plotted the convergence of the period 
and amplitude parameters of the fifth signal to this solution in 
Fig. |7] and the posterior probability densities of these two pa- 
rameters in Fig. [8] This signal was significantly present in the 
data and the five-Keplerian model was found to have a posterior 
probability of 460 times that of the four-Keplerian one, which 
exceeds the detection threshold of 150. We note that around the 
2000th member of the Markov chains in Fig. [7] (bottom panel), 
the variance of the proposal density has been decreased heavily 
to increase the convergence rate because the adaptation of the 
proposal density to the posterior is very low when the chain dis- 
covers a narrow probability maximum in the parameter space, as 
happens in this particular case because the width of the poste- 
rior density of the period parameter is much narrower than the 
corresponding parameter space (see Fig. [7] top panel). 

Samplings of both five- and six-Keplerian models revealed 
an additional probability maximum in the period space between 
160-200 days with two maxima at roughly 170 and 190 days, 
but the posterior samplings did not constrain the amplitude of 
this signal from below, which violates the detection criteria. 
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Fig. 7. Convergence of three Markov chains to a signal with a 
period of 43.7 days and an amplitude of 0.8 ms _1 with randomly 
selected initial states (arrows) in the period space between 12.9 
and 66.7 days given the combined HARPS and HIRES data set 
and a pure white noise model. 
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Fig. 8. Posterior probability densities of the period and ampli- 
tude of the fifth most significant (fourth shortest periodicity) sig- 
nal at a period of 43.7 days given the combined HARPS and 
HIRES data set and a pure white noise model. 



Therefore, using the white noise model leads to the conclusion 
that there are five periodic signals in the combined HARPS and 
HIRES data. The fact that the 190 day signal observed in the 
HARPS data alone did not satisfy the detection criteria is un- 
fortunate. While it did appear as a probability maximum given 
the combined data as well, it could not be constrained as it was 
for the HARPS data alone. This suggests that this signal is not a 
genuine Doppler signal of planetary origin but caused by either 
pure chance in the HARPS velocities or by some (pseudo) pe- 
riodic source of noise or bias in the HARPS data not accounted 
for by the pure white noise model. 

On the other hand, the nature of the 43.7 day signal in the 
combined data set is unknown and we cannot rule out its origin 
as a genuine Doppler signature of planetary origin. It is always 
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possible that a signal cannot be detected independently in any of 
the available datasets but only in the combined set because of the 
better phase-coverage, longer baseline, and a greater number of 
measurements in the combined set. Yet, it could also be caused 
by systematic differences between the different data sets, such 
as small instrumental biases that give rise to weak signals in the 
combined data set when not accounted for. However, we do not 
find any indication of such biases and in fact the HARPS and 
HIRES data sets cannot be shown to contradict with one another 
given the five-Kepleri an model bas e d on th e Bayesian model in- 
adequacy criterion of iTuomi et al.l d201 ll) . Further, it might be 
a coincidence that the fifth signal appears at a period that, if of 
planetary origin, implies a 3:2 commensurability of the corre- 
sponding planet with GJ 581 d that has an orbital period of 66.7 
days. Such an orbital configuration would likely enable the sta- 
bility of the system in long term but, while this is compatible 
with a hypothesis that the 43.7 day signal is caused by a fifth 
planet candidate in the system, it does not imply the existence 
of five planets around GJ 581 unless the existence of the fifth 
signal is supported by future data. Whether this explanation for 
the 43.7 day signal is valid or not also remains to be investigated 
with better noise models. 

We note that the excess noise levels (estimates of parameter 
cr y ) of the HARPS and HIRES velocities were 1.50 [1.15, 1.86] 
and 2.20 [1.63, 2.90] ms _1 , respectively, given the best model 
with k — 5. 

4.2. Red noise model 

We tested the dependence of the obser ved signals on the se- 
lected noise model and the conclusions of lBaluevl (1201 2l) by tak- 
ing into account red features in the radial velocity noise using 
a moving average mod el (M A). This model was constructed as 
in ITuomi et al.l (I2~012al) and [juomi et al. (2012b) but had only 
one component, making the model a first order MA model with 
Gaussian white excess noise. Mathematically, this noise model 
can be written as 

m = r k (ti) + y + €j + (p[mi-i - r k (ti-i) - y\ exp [a(t^ - ?,■)], (3) 

where each measurement m, at epoch f, is described using the su- 
perposition of k Keplerian signals (r r ); a reference velocity (y); 
Gaussian random variable corresponding to the white noise in 
the data (e,) with zero mean and variance cr? + cr y , where cr,- is 
the estimated instrument uncertainty from the spectral reduction 
pipeline and 07 is an unknown excess noise component some- 
times referred to as "stellar jitter"; parameter <p is a parameter 
describing the MA magnitude and the term in brackets corre- 
sponds to an exponential decrease of this magnitude as a func- 
tion of time difference of the ;th a nd i - 1th me asurements in 
accordance with the formulation o f lBaluevl(l20T2h . Parameter a 
describes the timescale of the correlations in the measurement 



As observed by iBaluevI d2012l) . the HIRES data alone ac- 
tually supports the existence of the planets GJ 581 b, c, and e 
despite the fact that only b and c could be detected usi ng a pure 
white noise model and Bayesian analyses of the data dGregorvi 
1201 ll) . Using the MA model described above, the samplings of 
the parameter space of the three-Keplerian model showed be- 
yond reasonable doubt t hat these three s ignals are present in the 
122 HIRES velocities of Vog t et all (120 1 Oh by increasing the pos- 
terior probability of the model by a factor of 8.2x10 greater 
than that of the two-Keplerian model that corresponded to the 
two signals at 5.4 and 12.9 days. The third signal with a period 



ode; 3.1490660 0.20828876 
; 3.1 490884 /a 4 : 0.736341 1 8 

4.04438761 E-04 



□de: 1.7272077 
: 1.6647037 
: 0.30832356 



P e [day] 



Fig. 9. As in Fig. [T]but for the velocity amplitude and period 
of GJ 581 e as obtained from the 122 HIRES velocities of 
IVogt et al.1 (l2010l) data using the MA model. 



of 3.1 days was constrained well and found consistent with the 
corresponding signal obtained from the HARPS velocities. We 
show the posterior densities of the period and velocity ampli- 
tude of GJ 581 e in Fig.|9]to demonstrate the significance of this 
detection. 

We could not spot a fourth significant periodic signal in the 
HIRES data using the posterior samplings and the MA model. 
However, we performed samplings anyway and spotted three 
probability maxima for the period of the fourth signal at roughly 
70, 300, and 500 days. However, we could not constrain any so- 
lutions for the four-Keplerian model and essentially ended up 
drawing samples from the prior of the period, though the chains 
were attracted to these three maxima throughout the samplings. 

Using the MA model also decreased the estimated amount of 
white noise in the data. While this white noise was found to be 
2.20 [1.63, 2.90] ms _1 as estimated with the parameter 07 in the 
pure white noise model, it decreased to 1.81 [1.37, 2.72] ms"' 
for the MA model. 

The HARPS velocities, however, did appear to favour a 
model with five periodic signals. We spotted a 190 day peri- 
odicity in the HARPS velocities using the MA model of Eq. 
(0. Compared to the four-Keplerian model, the five-Keplerian 
one had 540 times greater posterior probability and the pe- 
riod and amplitude of the signal at 190 days were well con- 
strained in accordance of our detection criteria. Therefore, the 
HARPS velocities, as was the case for the analyses with the 
pure white noise model, supports the existence of five indepen- 
dent sour ces of periodicity c orresponding to the four planets re- 
ported in Mayor et al. (2009) and one previously unknown signal 
at roughly 190 day period. We plotted the densities of parame- 
ters P and K of this signal in Fig. [10] (top panels) together with 
the corresponding phase-folded signal with the four planetary 
signatures of GJ 581 b, c, d, and e, and the MA component of 
the noise removed (Fig. [10] bottom panel). The white noise pa- 
rameter cry decreased from 1.50 [1.15, 1.86] to 1.19 [0.98, 1.64] 
ms"' when adding the MA component to the statistical model. 
This indicates, as was found for the HIRES data, that there are 
significant correlations in the noise of the HARPS data as well 
and they have to be accounted for when modelling the velocities. 

When analysing the combined HARPS and HIRES data set 
with the MA model of Eq. (0J, we did not have any trouble recov- 
ering the signals of the four planets orbiting the star. However, 
as was the case with the white noise model, we failed to iden- 
tify the existence of the 190 day signal as a significant one when 
sampling the period space above 67 days with a five-Keplerian 
model. We did observe a probability maximum at 190 days in 
the period space but could not constrain the corresponding sig- 
nal from above and below in the period space and below in 
the amplitude space in violation of the detection criteria. Also, 
we did not find a signal at a period of 43.7 days, which indi- 
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Fig. 10 . As in Fig.[3]but for the fifth signal in the lForveille et alJ 
d201 lb data given the MA model of Eq. (0. 



cates that the corresponding periodicity detected using a pure 
Gaussian white noise model was in fact most likely a spurious 
signal generated by insufficient noise modelling. This conclu- 
sion is supported by the fact that for k = 4 the MA model was 
found to have a Bayesian evidence of -797.0 on the logarithmic 
scale, whereas the corresponding Bayesian evidence of the white 
noise model was -813.1. Further, the log-Bayesian evidence of 
the white noise model with k — 5 was only -806.3, which im- 
plies that the MA model with k — 4 is the best description of 
the data out of the set of models tested. We show the parameter 
MAP estimates and their 99% BCSs of this best model in Table 

m 

We note that samplings of a six-Keplerian model did not 
show any evidence in favour of additional periodic signals in 
the combined data set. 



5. Discussion 

Estimating the number of planetary companions orbiting GJ 
581, or equivalently, estimating the parameter k favoured by the 
data for this star, has b een attempted by several authors over 
the past few years (e.g. iBonfils et all 120051: lUdrv et all [2 007: 



May or et all 120091: IVogt et all 1201 ft iTuomil 1201 U iGregorv , 
l201ltrVogt et all[2012l : lBaluevll2012h . Improvements in statisti- 
cal analysis techniques have proven efficient means of improving 
this estimate, i.e. especially showing that some signals have actu- 
ally been false positives, but the most significant developments 
have relied on the availability of new high-precision data. We 
have assessed the number of significant periodic signals in the 
combined HARPS and HIRES data of GJ 581 and conclude that 
there is evidence in favour of four signal s, corresponding to the 
four planets orbiting the star observed by Mavor et al] d2009l) . 

By re-analysing sequentially the HARPS data after each 
observing period, i.e. each full year, we have shown that the 
Bayesian statistical techniques were at least as sensitive in de- 
tecting periodic signals in the data as the more commonly ap- 
plied periodogram-based analysis tools (TableQ]). Moreover, we 
did not find any significant periodicities in these HARPS data 
subsets that did not show as significant signals when the amount 
of data was increased. This implies that the Bayesian detection 



criteria (Tuomil l2012h are not sensitive to false positives aris- 
ing from noise by chance in practice and any periodic signals 
that satisfy these criteria are likely genuine signals present in the 
data. 

However, insufficient modelling and possible biases in any 
single data set can produce false positives that might be inter- 
preted as planetary signals and could contaminate the statis- 
tics of exoplanet properties. For instance, the HARPS data was 
found to favour five periodicities after the sixth full observing 
period and 197 velocities. This fifth signal at a period of 190 
days did not disappear after the inclusion of another 43 veloc- 
ity measurements a nd wa s clearly present in the data published 
by iForveille et al.l d201 ll) . However, another independent data 
set obtained using the HIRES did not support the existence of 
this signal. This was evident when the signal could not be de- 
tected significantly in the combined HARPS and HIRES data 
set. Therefore, we conclude that the 190 day periodicity in the 
HARPS data set, while being significantly present, was likely 
not of planetary origin but a spurious periodicity caused by some 
source of (systematic) noise that had variations resembling peri- 
odic behaviour. Alternatively, it is possible that the noise model 
does not describe the HIRES and/or HARPS data well enough, 
which might result in the inability to detect this signal in the 
combined set. However, we do not see any evidence of such in- 
adequate modelling as the two data sets are consi stent accord- 
ing to the Bayesian model inadequacy criterion of iTuomi et al.l 
(1201 lh . 

The combined HARPS and HIRES data set, however, did 
contain a fifth periodic signal. This time the weak fifth signal 
occurred at a period of 43.7 days. This result was obtained by 
modelling the two velocity sets using Gaussian white noise mod- 
els that are almost a standard way of describin g radia l velocities 
in practice. As was shown recently by Baluevl d2012l) . the radial 
velocity noise of GJ 581 is not white but contains significant 
correlations (red noise) that, if not accounted for, could mimic 
periodic behaviour. We found this to be the case and observed 
that the statistical model containing a first order moving average 
component improved the model significantly and enabled us to 
describe the combined data the best with only four periodic sig- 
nals. Therefore, the 43.7 day signal observed in the combined 
data set using the white noise model is the most likely a spurious 
signal caused by insufficient modelling and consequent misin- 
terpretation of the data and is not a genuine Doppler signal of 
planetary origin. This does not, however, imply beyond reason- 
able doubt that the 43.7 day periodicity is not a genuine signal 
because it is possible that the noise model is still not optima^ 
and that a better model could enable the recovery of the signal. 

These results show that not only can insufficient modelling 
lead to detections of false positives, but relying on a single data 
set can also lead to detections of signals that are not supported by 
other independent measurements. False positives can be caused 
by these two sources in practice and care is needed to assess 
whether any given weak signal in any single data set is caused 
by either of these two effects or not. However, our results sug- 
gest that while the Bayesian detection criteria are not particu- 
larly sensitive to false positives, they can fall a victim of these 
two sources of false positives. 

Finally, based on our results, it looks like there are four and 
only four planets orbitin g GJ 581. This means that the 32 day 
signal observed by Vogt et al.l d2012l) was a false positive, though 



3 Actually this is not only possible but rather certain because the ori- 
gin, nature, time-evolution, and other features of radial velocity noise 
are very poorly understood at the moment. 
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Table 2. Four-Keplerian solution of the HARPS and HIRES velocities of GJ 581 together with the parameters quantifying the noise 
properties. MAP estimates of the parameters and their 99% BCSs. 



Parameter 


e 


b 


c 


d 


P [days] 


3.14912 [3.14865, 3.14963] 


5.36863 [5.36840, 5.36884] 


12.9169 [12.9104, 12.9235] 


66.63 [66.31,66.92] 


e 


0.16 [0, 0.36] 


0.03 [0, 0.05] 


0.03 [0, 0.20] 


0.17 [0, 0.51] 


K [ms -1 ] 


1.72 [1.36, 2.11] 


12.63 [12.24, 13.05] 


3.21 [2.64, 3.73] 


1.87 [1.12, 2.47] 


a> [rad] 


0.6 [0, 2/r] 


0.9 [0, 2n\ 


1.3 [0, 2/r] 


6.2 [0, 2n\ 


M [rad] 


5.8 [0, 2tt] 


5.9 [0, 2n] 


1.9 [0, 2tt] 


0.6 [0, 2n] 


a [AU] 


0.0285 [0.0266, 0.0302] 


0.0407 [0.0380, 0.0431] 


0.0731 [0.0684, 0.0773] 


0.218 [0.204, 0.231] 


m p sin i [M ffi ] 


1.8 [1.3,2.3] 


15.9 [13.7, 17.9] 


5.4 [4.2, 6.5] 


5.2 [3.3,7.5] 


Parameter 


HARPS 


HIRES 






y [ms~'] 


0.11 [-0.46, 0.61] 


0.14 [-0.87, 1.22] 






cry [ms -1 ] 


1.23 [1.06, 1.72] 


1.62 [1.15, 2.32] 






a [day -1 ] 


0.057 [0, 0.446] 


0.120 [0, 0.863] 






<P 


0.49 [0.26, 0.95] 


0.91 [0.45, 1] 







the evidence in favour of this signal in their analyses w as mod- 
erate a t most anyway. Also, the periodogram analyses of lBaluevl 
(120121) obtained a decreased significance for the companion GJ 
58 1 d because of the inconsistensie s in the periodo gram analyses 
described in Section lZSland in e.g. | Tupmil d2012l) . However, the 
HARPS velocities of F orveille et al.1 (1201 11) obtained from the 
spectra us ing the cross -correlation function method are not opti- 
mal dPepe et al.Ll2002l) but their accuracy can be significantly im- 
proved when using more sophisticat ed spectral reduction meth- 
ods, such as the TERRA algorithm dA nglada-E scude & Butlerl 
l2012t lAnglada -Escude & Tuoml 120121) . Improving the noise 
modelling can also lead to more accurate descriptions of the ve- 
locities and therefore enable the detections of one or some of the 
signals that were only obtained from the HARPS data or with a 
particular noise model at 190 and 44 day periods. A better un- 
derstanding of the noise inherent i n measuring RVs fo r M-dwarf 
stars is essential because curre nt (Barn eset al.L I2012D and future 
( Ramsev et all I2008L iMahadevan et al.Ll2009t) planet search sur- 
veys aim to focus their attenti on on the inactive a nd slowly rotat- 
ing subset of these stars (see lJenkins et al.Ll2009l) . For these rea- 
sons, and because GJ 581 is still a target of e.g. the HARP S and 
HIRES exoplanet surveys (according to lVogt et alJ d2012l) there 
is already a considerable amount of additional HIRES data), we 
do not expect that the debate over the number of k for GJ 581 
is over. This is especially interesting because the period-space 
between GJ 581 c and d is known to contain dynamically sta- 
ble orbits for low-mass planets and happens to correspond to the 
liquid water habitable zone of the star. 
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